Ramalho et al 2021

Usefull librarys.

library(rmarkdown)    # You need this library to run this template.
library(epuRate)      # Install with devtools: install_github("holtzy/epuRate", force=TRUE)
library(DT)           # pour faire de joies affichage de dataframe
library(readtext)
#library(tcltk)   # choose directory
library(readxl)
library(openxlsx)
library(ggplot2)
library(tidyverse)

# detach("package:MASS")

library(afex)       # nouveau package, surveiller les MAJ
library(lme4)
library(RVAideMemoire)
library(LMERConvenienceFunctions)
library(phia)
library(multcomp)
library(lmerTest)
library(emmeans)
library(ggpubr)
library(ggeffects)
library(splines)
library(DHARMa)
library(glmmTMB)
library(sjPlot) 
library(sjmisc) 

library(rstanarm)
library(brms)  # for models
library(bayesplot)
library(bayestestR)

Ramalho B. L., Moly J., Raffin E., Bouet R., Harquel S., Farnè A., Reilly K.T. (2021) Electrocutaneous stimulation to the face inhibits motor evoked potentials in the hand: Face-hand sensorimotor interactions revealed by afferent inhibition. EJN

https://pubmed.ncbi.nlm.nih.gov/34796553/

In this study we used the afferent inhibition protocol to assess sensorimotor interactions between the hand and body parts represented close to the hand area in the sensorimotor cortex: the face and the upper limb. We found that the amplitude of MEPs in the right FDI was reduced when the TMS pulse was preceded by an electrocutaneous stimulus on the face (upper lip or cheek) or on the upper limb (arm or forearm), but this inhibition was more robust when associated with face stimulation. These results provide the first evidence for afferent inhibition between the skin on the face and a hand muscle, but also between the skin on the upper limb and a hand muscle.

Key words:
Time series; GLM; Gamma; Splines; Power; excel sheets;

Datas

Extract

Here we define a function to extract datas from exel files:
- Read excel file
- Explore all sheet (coresponding to subjects)
- TP stimuation become O ms timing.
Like this I can compare each time point (15,15,35….) with 0 and we keep TP variance.
- Each mesure was centered (\(Mesure.TP = Mesure / TP\))
Exitability difference with Test Pulse baseline.

We add one factor : Localization
LIP & CHEEK : loc1
FOREAM & ARM : loc2

read.all_subject <- function(prefix){
  
  filename <- paste("/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Datas_raw/", prefix, "_all_subjects.xlsx", sep = "")
 
  sheets.names <- excel_sheets(filename)
  datas.sheets <- lapply(excel_sheets(filename), read_excel, path = filename)
  
  datas.raw <- data.frame()
  # first sheet is MEAN
  for (xi_suj in 1:length(datas.sheets)){
      datas.raw <- bind_rows(datas.raw, mutate(datas.sheets[[xi_suj]], Subject = sheets.names[xi_suj]) )
  }
  rm(xi_suj)
  rm(filename, sheets.names, datas.sheets)
  
  # Extract subject TP mean 
  datas.raw %>%
    dplyr::select(Subject, TP) %>%
    group_by(Subject) %>%
    summarise(TP_mean = mean(TP),
              TP_sd = sd(TP),
              TP_med = median(TP)) -> datas.TP
  
  
  datas.raw %>%
    mutate("0" = TP) %>%
    dplyr::select(-TP) %>%
    gather(key = "Time", value = "Mesure", -Subject) %>%
    mutate(Time = as.factor(Time)) %>%
    drop_na() %>%
    left_join(datas.TP) %>%
    group_by(Subject, Time, Mesure) %>%
    data.frame() %>%
    mutate(Mesure.TP = Mesure/TP_mean,
           Mesure.TP.med = Mesure/TP_med) %>%
    mutate(Body = as.factor(prefix)) %>%
    mutate(Subject = as.factor(Subject)) -> datas.raw
  
  return(datas.raw)
}

datas.forearm <- read.all_subject("FOREARM")
datas.arm <- read.all_subject("ARM")
datas.cheek <- read.all_subject("CHEEK")
datas.lip <- read.all_subject("LIP")

bind_rows(datas.forearm, datas.arm, datas.cheek, datas.lip) %>%
  mutate(Subject = as.factor(Subject),
         Body = as.factor(Body)) %>%
  dplyr::select(Body, Subject, Time, Mesure, Mesure.TP, Mesure.TP.med) -> datas.raw

# Add localization factor
datas.raw %>%
  mutate(Localization = case_when(Body == "LIP" ~ "Loc1",
                                  Body == "CHEEK" ~ "Loc1",
                                  Body == "FOREARM" ~ "Loc2",
                                  Body == "ARM" ~ "Loc2")) %>%
  mutate(Localization = as.factor(Localization)) -> datas.raw



datas.raw %>%
  dplyr::select(-Subject) %>%
  gtsummary::tbl_summary(by = "Body") %>%
  gtsummary::add_p()
Characteristic FOREARM, N = 2,0141 ARM, N = 1,7201 CHEEK, N = 1,5181 LIP, N = 1,8241 p-value2
Time >0.9
0 476 (24%) 410 (24%) 346 (23%) 416 (23%)
15 195 (9.7%) 166 (9.7%) 147 (9.7%) 178 (9.8%)
25 190 (9.4%) 166 (9.7%) 148 (9.7%) 179 (9.8%)
35 185 (9.2%) 163 (9.5%) 152 (10%) 180 (9.9%)
45 193 (9.6%) 159 (9.2%) 143 (9.4%) 172 (9.4%)
55 193 (9.6%) 165 (9.6%) 144 (9.5%) 167 (9.2%)
65 190 (9.4%) 163 (9.5%) 147 (9.7%) 172 (9.4%)
75 199 (9.9%) 165 (9.6%) 142 (9.4%) 180 (9.9%)
85 193 (9.6%) 163 (9.5%) 149 (9.8%) 180 (9.9%)
Mesure 0.78 (0.39, 1.37) 0.87 (0.46, 1.36) 0.58 (0.24, 1.17) 0.58 (0.27, 1.13) <0.001
Mesure.TP 0.81 (0.41, 1.23) 0.88 (0.46, 1.31) 0.72 (0.26, 1.36) 0.72 (0.35, 1.35) <0.001
Mesure.TP.med 0.90 (0.49, 1.38) 0.90 (0.48, 1.33) 0.83 (0.32, 1.55) 0.94 (0.51, 1.75) <0.001
Localization <0.001
Loc1 0 (0%) 0 (0%) 1,518 (100%) 1,824 (100%)
Loc2 2,014 (100%) 1,720 (100%) 0 (0%) 0 (0%)
1 n (%); Median (IQR)
2 Pearson's Chi-squared test; Kruskal-Wallis rank sum test 
rm(datas.arm, datas.cheek, datas.forearm, datas.lip, read.all_subject)

Display

Distribution

Some displays to explore distributions. Here time is categorical factor.

datas.raw %>%
  mutate(Time = Time) %>%
  ggplot(aes(x = Time, y = Mesure.TP, color = Time)) +
    geom_jitter(color = "black", alpha = 0.1) + 
    geom_violin(alpha = 0, position = position_dodge(width = 0.9)) + 
    geom_boxplot(width=.2, alpha = 0, outlier.alpha = 0, 
              position = position_dodge(width = 0.9)) +
    stat_summary(fun.y = "mean", 
               geom = "point", 
               size = 2, 
               position = position_dodge(width=0.9), 
                               color = "black") + 
    coord_cartesian(ylim =c(-0.1, 5)) +
    facet_wrap(~Body)

Datas does’t looks Gaussian.

Focus on “ARM” :

datas.raw %>%
  mutate(Time = Time) %>%
  filter(Body == "ARM") %>%
  ggplot(aes(x = Time, y = Mesure.TP, color = Time)) +
    geom_violin(alpha = 0, position = position_dodge(width = 0.9)) + 
    geom_boxplot(width=.2, alpha = 0, outlier.alpha = 0, 
              position = position_dodge(width = 0.9)) +
    stat_summary(fun.y = "mean", 
               geom = "point", 
               size = 2, 
               position = position_dodge(width=0.9), 
                               color = "black") + 
    coord_cartesian(ylim =c(-0.1, 4)) +
    facet_wrap(~Subject)

As we can see. It’s usefull to add a “0” time as a TP periode. Variance within TP condition help to see a difference between another time.

Variability

Explore inter-subject variability

datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) + 
    geom_errorbar(aes(ymin = Mesure.TP-ci, ymax = Mesure.TP+ci), 
                  colour = "black", alpha = 0.4, 
                  width = .1, position = pd) +
    geom_line(position = pd) +
    geom_point(position = pd, size = 3) +
    ggtitle("Subject average")

datas.raw %>%
    group_by(Subject, Body, Time) %>%
    summarise(Mesure.TP = mean(Mesure.TP)) %>%
    mutate(Time = as.numeric(as.character(Time))) %>%
    ggplot(aes(x = Time, y = Mesure.TP, color = Subject)) +
      geom_line(width=.2, show.legend = FALSE) +
    facet_wrap(~Body, scales = "free") +
    theme(legend.position="bottom")+
    ggtitle("Subject datas")

High Subject’s heterogeneity for Mesure.TP along Time.

Stats

I did 2 analysis:

  • like the Karen’s paper, I use Time as a categorial variable
  • To introduce a inter-dependance between time periode, also I use Time as a numerical variable

We use a glm method because we have to:
- Consider the Gamma distribution (Datas are continuous and bounded) - We accound of inter-Subject variability and Time by Body variability (random effect).

Time as factor

First, Time is consider as factor.
So we explore 2 fixe effects and there interation:
- Body (4 levels)
- Time (9 levels)

We consider 3 random effets:
- Subject intercept
- Body intercept - Time slope by Body

Models

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Fact.RData")
# model.03 <- glmer(Mesure.TP ~ (Body + Time)^2 + (1|Subject) + (Time|Body),
#                   family = Gamma("identity"),
#                   data = datas.raw)

# model.03 <- glmer(Mesure.TP ~ (Body + Time)^2 + (1|Subject) + (Time|Body) + (1|Localization),
#                   family = Gamma("identity"),
#                   data = datas.raw)

#save(model.01, model.02, model.03,
#      file = "/Volumes/crnldata/dycog/Epilepto/Karen/TMS/Scripts/Models_Time_Fact.RData")

Anova(model.03, type = "II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Mesure.TP
##             Chisq Df Pr(>Chisq)    
## Body       10.902  3   0.012270 *  
## Time      132.020  8  < 2.2e-16 ***
## Body:Time  57.385 24   0.000148 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We have a significative interaction Body:Time.

Post-hoc

Effect’s displays for raw datas and estimate mesures.

datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) + 
    geom_errorbar(aes(ymin = Mesure.TP-ci, ymax = Mesure.TP+ci), 
                  colour = "black", alpha = 0.4, 
                  width = .1, position = pd) +
    geom_line(position = pd) +
    geom_point(position = pd, size = 3) +
    ggtitle("Raw datas")

plot_model(model.03,
           type = "pred",
           terms = c("Time", "Body")) +
    ggtitle("Estimate datas")

Here time-pairwise by Body, Bonferroni correction:

em <- emmeans(model.03, ~ Time | Body)
#CLD(em)
# contrast(em, interaction = c( "pairwise"), adjust = "bonferroni")
des_cont <- list("0_15" = c(-1, 1, 0, 0, 0, 0, 0, 0, 0),
                 "0_25" = c(-1, 0, 1, 0, 0, 0, 0, 0, 0),
                 "0_35" = c(-1, 0, 0, 1, 0, 0, 0, 0, 0),
                 "0_45" = c(-1, 0, 0, 0, 1, 0, 0, 0, 0),
                 "0_55" = c(-1, 0, 0, 0, 0, 1, 0, 0, 0),
                 "0_65" = c(-1, 0, 0, 0, 0, 0, 1, 0, 0),
                 "0_75" = c(-1, 0, 0, 0, 0, 0, 0, 1, 0),
                 "0_85" = c(-1, 0, 0, 0, 0, 0, 0, 0, 1))
contrast(em, des_cont, adjust = "bonferroni")
## Body = FOREARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15       0.0454 0.0786 Inf   0.577  1.0000
##  0_25      -0.0324 0.0745 Inf  -0.435  1.0000
##  0_35      -0.1218 0.0694 Inf  -1.754  0.6359
##  0_45      -0.1484 0.0686 Inf  -2.164  0.2435
##  0_55      -0.1356 0.0692 Inf  -1.959  0.4012
##  0_65      -0.1845 0.0668 Inf  -2.761  0.0460
##  0_75      -0.0254 0.0747 Inf  -0.340  1.0000
##  0_85      -0.0760 0.0720 Inf  -1.055  1.0000
## 
## Body = ARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15       0.1280 0.0861 Inf   1.486  1.0000
##  0_25      -0.0933 0.0738 Inf  -1.264  1.0000
##  0_35      -0.1907 0.0706 Inf  -2.700  0.0555
##  0_45      -0.3533 0.0621 Inf  -5.686  <.0001
##  0_55      -0.4000 0.0636 Inf  -6.293  <.0001
##  0_65      -0.2764 0.0715 Inf  -3.863  0.0009
##  0_75      -0.1191 0.0786 Inf  -1.516  1.0000
##  0_85      -0.0596 0.0824 Inf  -0.723  1.0000
## 
## Body = CHEEK:
##  contrast estimate     SE  df z.ratio p.value
##  0_15      -0.0981 0.0651 Inf  -1.507  1.0000
##  0_25      -0.1923 0.0628 Inf  -3.063  0.0175
##  0_35      -0.1439 0.0650 Inf  -2.214  0.2145
##  0_45      -0.0639 0.0690 Inf  -0.926  1.0000
##  0_55      -0.1945 0.0617 Inf  -3.151  0.0130
##  0_65      -0.1228 0.0653 Inf  -1.881  0.4803
##  0_75      -0.1221 0.0655 Inf  -1.865  0.4972
##  0_85      -0.1843 0.0630 Inf  -2.927  0.0274
## 
## Body = LIP:
##  contrast estimate     SE  df z.ratio p.value
##  0_15       0.0164 0.0736 Inf   0.222  1.0000
##  0_25      -0.0379 0.0724 Inf  -0.523  1.0000
##  0_35      -0.0592 0.0723 Inf  -0.819  1.0000
##  0_45      -0.2962 0.0603 Inf  -4.913  <.0001
##  0_55      -0.3599 0.0612 Inf  -5.885  <.0001
##  0_65      -0.3320 0.0593 Inf  -5.596  <.0001
##  0_75      -0.1459 0.0660 Inf  -2.210  0.2167
##  0_85      -0.0735 0.0700 Inf  -1.051  1.0000
## 
## P value adjustment: bonferroni method for 8 tests

First analysis found significant :
- LIP inhibition : 45, 55 and 65 ms
- CHEEK inhibition : 25, 55 and 85 ms
- FOREARM 65 ms
- ARM 45, 55 and 65 ms

Time as numeric

Is intuitive to consider a Time as numeric values, not factor. There is a dependence between levels.
Also, Time is not a linear effect. So we fit a non-linear Time effect with spline with 7 breakpoints.

Models

We use exactly the same model as previously.

datas.raw %>%
  mutate(Time = as.numeric(as.character(Time))) -> datas.raw.num

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume.RData")

# model.num.04 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),
#                        family = Gamma("identity"),
#                        data = datas.raw.num,
#                        method = "LRT")
 

 # save(model.num.04,
 #      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume.RData")
 

Anova(model.num.04$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##                     Chisq Df Pr(>Chisq)    
## (Intercept)        8.4606  1   0.003629 ** 
## Body               1.5617  3   0.668099    
## ns(Time, 7)      122.1927  6  < 2.2e-16 ***
## Body:ns(Time, 7)  54.1249 18  1.755e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model.num.04)
## Mixed Model Anova Table (Type 3 tests, LRT-method)
## 
## Model: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time, 
## Model:     7) | Body)
## Data: datas.raw.num
## Df full model: 66
##                  Df   Chisq Chi Df Pr(>Chisq)  
## Body             66  0.0052      0             
## ns(Time, 7)      60 16.5034      6    0.01129 *
## Body:ns(Time, 7) 48 19.2304     18    0.37776  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We always find a significant interaction Body:Time

Post-hoc

Effect’s displays for raw datas and estimate mesures.

datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) + 
    geom_smooth(width=.2, show.legend = TRUE, se = FALSE) +
    ggtitle("Raw datas")

plot_model(model.num.04$full_model,
           type = "pred",
           terms = c("Time", "Body")) +
    ggtitle("Estimate datas")

There is an effect, ok, but where ? In this case, When ?
To explore this issue, we estimate Mesure.TP for differents Time periode ((0, 15, 25, 35, 45, 55, 65, 75, 85). In this way, as a factor analysis, we compare each 8 time periode to 0 (baseline). Like this, we could express when Mesure.TP is significantly different from baseline.

em <- emmeans(model.num.04, ~ Time | Body, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)),
              type = "pred")


des_cont <- list("0_15" = c(1, -1, 0, 0, 0, 0, 0, 0, 0),
                 "0_25" = c(1, 0, -1, 0, 0, 0, 0, 0, 0),
                 "0_35" = c(1, 0, 0, -1, 0, 0, 0, 0, 0),
                 "0_45" = c(1, 0, 0, 0, -1, 0, 0, 0, 0),
                 "0_55" = c(1, 0, 0, 0, 0, -1, 0, 0, 0),
                 "0_65" = c(1, 0, 0, 0, 0, 0, -1, 0, 0),
                 "0_75" = c(1, 0, 0, 0, 0, 0, 0, -1, 0),
                 "0_85" = c(1, 0, 0, 0, 0, 0, 0, 0, -1))
contrast(em, des_cont, adjust = "bonferroni")
## Body = FOREARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15      0.11936 0.0592 Inf   2.016  0.3507
##  0_25      0.16543 0.0543 Inf   3.049  0.0184
##  0_35      0.15410 0.0569 Inf   2.707  0.0543
##  0_45      0.10241 0.0599 Inf   1.710  0.6982
##  0_55      0.14039 0.0547 Inf   2.568  0.0819
##  0_65      0.16410 0.0571 Inf   2.872  0.0326
##  0_75      0.10402 0.0663 Inf   1.569  0.9340
##  0_85      0.18853 0.0628 Inf   3.004  0.0213
## 
## Body = ARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15     -0.05211 0.0721 Inf  -0.723  1.0000
##  0_25      0.04108 0.0603 Inf   0.681  1.0000
##  0_35      0.12770 0.0621 Inf   2.058  0.3171
##  0_45      0.12296 0.0644 Inf   1.910  0.4493
##  0_55      0.16844 0.0572 Inf   2.943  0.0260
##  0_65      0.16148 0.0630 Inf   2.563  0.0830
##  0_75      0.03861 0.0723 Inf   0.534  1.0000
##  0_85      0.07519 0.0722 Inf   1.041  1.0000
## 
## Body = CHEEK:
##  contrast estimate     SE  df z.ratio p.value
##  0_15     -0.08051 0.0739 Inf  -1.089  1.0000
##  0_25      0.03605 0.0635 Inf   0.568  1.0000
##  0_35      0.21975 0.0626 Inf   3.508  0.0036
##  0_45      0.35782 0.0588 Inf   6.081  <.0001
##  0_55      0.38872 0.0571 Inf   6.806  <.0001
##  0_65      0.29071 0.0648 Inf   4.483  0.0001
##  0_75      0.11239 0.0778 Inf   1.444  1.0000
##  0_85      0.05736 0.0829 Inf   0.692  1.0000
## 
## Body = LIP:
##  contrast estimate     SE  df z.ratio p.value
##  0_15      0.00231 0.0667 Inf   0.035  1.0000
##  0_25      0.00442 0.0595 Inf   0.074  1.0000
##  0_35      0.08749 0.0624 Inf   1.402  1.0000
##  0_45      0.27856 0.0567 Inf   4.912  <.0001
##  0_55      0.37798 0.0523 Inf   7.223  <.0001
##  0_65      0.32072 0.0560 Inf   5.727  <.0001
##  0_75      0.15236 0.0643 Inf   2.368  0.1430
##  0_85      0.07082 0.0703 Inf   1.007  1.0000
## 
## P value adjustment: bonferroni method for 8 tests

For instance, here we see a significative decrease for LIP at 45, 55 and 65 (p.value was control for multiple comparaison by bonferroni method).

Supplementary analysis

Following the advice of the expert reviewers from this first round of reviews we made substantial changes. All these analysis was recommended by reviewers.

Experiment Factor

We keep model with Time as numeric.
We add factor to code the experiment factor.

LIP & CHEEK : Experiment 1
FOREAM & ARM : Experiment 2

Because we are interst by fit the homogeneity within this Experiment factor.
We use this factor as random intercept.

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Experiment.RData")

# model.num.05 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 + 
#                        (1|Subject) + (ns(Time,7)|Body) + (1|Localization), 
#                        family = Gamma("identity"),
#                        data = datas.raw.num, 
#                        method = "LRT")
# 
# save(model.num.05,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Experiment.RData")

First we explore the difference between previous model.

anova(model.num.04, model.num.05)
## Data: data
## Models:
## model.num.04: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time, 
## model.num.04:     7) | Body)
## model.num.05: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time, 
## model.num.05:     7) | Body) + (1 | Localization)
##              npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)
## model.num.04   66 12683 13136 -6275.4    12551                     
## model.num.05   67 12685 13145 -6275.4    12551 0.0029  1     0.9572

There is no significative differences between these models. And we see an BIC and AIC increase with the last model. So we conclude to the lack of relevance for this new model.

Without surpise, we found the same significativity:

Anova(model.num.05$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##                     Chisq Df Pr(>Chisq)    
## (Intercept)        8.4759  1   0.003599 ** 
## Body               1.5386  3   0.673392    
## ns(Time, 7)      122.2503  6  < 2.2e-16 ***
## Body:ns(Time, 7)  54.0876 18  1.778e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.05$full_model, terms = c("Time", "Body")))

em <- emmeans(model.num.05$full_model, ~ Time | Body, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

contrast(em, des_cont, adjust = "bonferroni")
## Body = FOREARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15      0.11904 0.0592 Inf   2.010  0.3551
##  0_25      0.16529 0.0542 Inf   3.048  0.0184
##  0_35      0.15415 0.0569 Inf   2.710  0.0539
##  0_45      0.10234 0.0599 Inf   1.709  0.6991
##  0_55      0.14000 0.0547 Inf   2.561  0.0835
##  0_65      0.16375 0.0571 Inf   2.867  0.0332
##  0_75      0.10408 0.0663 Inf   1.570  0.9305
##  0_85      0.18800 0.0627 Inf   3.000  0.0216
## 
## Body = ARM:
##  contrast estimate     SE  df z.ratio p.value
##  0_15     -0.05128 0.0721 Inf  -0.711  1.0000
##  0_25      0.04127 0.0603 Inf   0.684  1.0000
##  0_35      0.12735 0.0621 Inf   2.051  0.3220
##  0_45      0.12308 0.0644 Inf   1.911  0.4478
##  0_55      0.16879 0.0572 Inf   2.949  0.0255
##  0_65      0.16176 0.0630 Inf   2.568  0.0819
##  0_75      0.03888 0.0722 Inf   0.538  1.0000
##  0_85      0.07571 0.0720 Inf   1.052  1.0000
## 
## Body = CHEEK:
##  contrast estimate     SE  df z.ratio p.value
##  0_15     -0.08016 0.0739 Inf  -1.085  1.0000
##  0_25      0.03640 0.0635 Inf   0.573  1.0000
##  0_35      0.21999 0.0626 Inf   3.513  0.0035
##  0_45      0.35803 0.0588 Inf   6.085  <.0001
##  0_55      0.38905 0.0571 Inf   6.814  <.0001
##  0_65      0.29121 0.0648 Inf   4.493  0.0001
##  0_75      0.11289 0.0778 Inf   1.451  1.0000
##  0_85      0.05732 0.0826 Inf   0.694  1.0000
## 
## Body = LIP:
##  contrast estimate     SE  df z.ratio p.value
##  0_15      0.00204 0.0667 Inf   0.031  1.0000
##  0_25      0.00451 0.0595 Inf   0.076  1.0000
##  0_35      0.08774 0.0624 Inf   1.407  1.0000
##  0_45      0.27841 0.0567 Inf   4.909  <.0001
##  0_55      0.37806 0.0523 Inf   7.224  <.0001
##  0_65      0.32106 0.0560 Inf   5.734  <.0001
##  0_75      0.15246 0.0643 Inf   2.370  0.1424
##  0_85      0.07143 0.0702 Inf   1.017  1.0000
## 
## P value adjustment: bonferroni method for 8 tests

Two separate analysis

Reviewers ask to analyse Experiment 1 and Experiment 2 datas separatly. They are affraid that our GLMM hide an Experiment effect.
So, we split (physically) datas according to Experiment 1 and Experiment 2 and we fit 2 models:

Experiment 1 (LIP & CHEEK)

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Supp.RData")

# model.num.06 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 + 
#                        (1|Subject) + (ns(Time,7)|Body), 
#                        family = Gamma("identity"),
#                        data = filter(datas.raw.num, Localization == "Loc1"), 
#                        method = "LRT")

Anova(model.num.06$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##                    Chisq Df Pr(>Chisq)    
## (Intercept)       1.5181  1     0.2179    
## Body              0.1043  1     0.7467    
## ns(Time, 7)      69.3255  6  5.622e-13 ***
## Body:ns(Time, 7)  3.8732  6     0.6938    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.06$full_model, terms = c("Time", "Body")))

em <- emmeans(model.num.06$full_model, "Time", 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15      -0.0333 0.0535 Inf  -0.622  1.0000
##  0_25       0.0289 0.0466 Inf   0.620  1.0000
##  0_35       0.1628 0.0473 Inf   3.446  0.0046
##  0_45       0.3238 0.0441 Inf   7.349  <.0001
##  0_55       0.3811 0.0414 Inf   9.201  <.0001
##  0_65       0.3015 0.0461 Inf   6.536  <.0001
##  0_75       0.1377 0.0538 Inf   2.559  0.0839
##  0_85       0.0747 0.0580 Inf   1.287  1.0000
## 
## Results are averaged over the levels of: Body 
## P value adjustment: bonferroni method for 8 tests

The same results as before

Experiment 2 (FOREAM & ARM)

# model.num.07 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 + 
#                        (1|Subject) + (ns(Time,7)|Body), 
#                        family = Gamma("identity"),
#                        data = filter(datas.raw.num, Localization == "Loc2"), 
#                        method = "LRT")

Anova(model.num.07$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##                      Chisq Df Pr(>Chisq)    
## (Intercept)         29.715  1  5.005e-08 ***
## Body              3054.248  1  < 2.2e-16 ***
## ns(Time, 7)      93350.704  6  < 2.2e-16 ***
## Body:ns(Time, 7) 23669.669  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.07$full_model, terms = c("Time", "Body")))

em <- emmeans(model.num.07$full_model, ~ Time, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15       0.0291 0.0338 Inf   0.859  1.0000
##  0_25       0.1284 0.0215 Inf   5.962  <.0001
##  0_35       0.1383 0.0137 Inf  10.099  <.0001
##  0_45       0.1204 0.0131 Inf   9.194  <.0001
##  0_55       0.1602 0.0135 Inf  11.861  <.0001
##  0_65       0.1601 0.0135 Inf  11.881  <.0001
##  0_75       0.0749 0.0132 Inf   5.695  <.0001
##  0_85       0.1357 0.0111 Inf  12.262  <.0001
## 
## Results are averaged over the levels of: Body 
## P value adjustment: bonferroni method for 8 tests
# save(model.num.06, model.num.07,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Supp.RData")

The same results as before exept for :
ARM 0_65 ms
FOREARM 0_55
becomes significant

Four separate analysis

The same as previously. Reviwers ask to analyse each body parts independently.
So we split datas according to Body levels and we fit 4 models:

LIP

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Four_separate.RData")


# model.num.08 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject), 
#                       family = Gamma("identity"), 
#                       data = filter(datas.raw.num, Body == "LIP"),
#                       method = "LRT")

Anova(model.num.08$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##               Chisq Df Pr(>Chisq)    
## (Intercept)  3.9684  1    0.04636 *  
## ns(Time, 7) 69.1470  6  6.116e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.08$full_model, terms = c("Time")))

em <- emmeans(model.num.08$full_model, ~ Time, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15      0.00165 0.0675 Inf   0.024  1.0000
##  0_25      0.00848 0.0598 Inf   0.142  1.0000
##  0_35      0.09688 0.0627 Inf   1.545  0.9795
##  0_45      0.28618 0.0579 Inf   4.939  <.0001
##  0_55      0.36836 0.0529 Inf   6.966  <.0001
##  0_65      0.30150 0.0578 Inf   5.213  <.0001
##  0_75      0.14755 0.0645 Inf   2.288  0.1772
##  0_85      0.07338 0.0709 Inf   1.035  1.0000
## 
## P value adjustment: bonferroni method for 8 tests

The same results as before

CHEEK

# model.num.09 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject), 
#                       family = Gamma("identity"), 
#                       data = filter(datas.raw.num, Body == "CHEEK"),
#                       method = "LRT")

Anova(model.num.09$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##              Chisq Df Pr(>Chisq)    
## (Intercept)  1.130  1     0.2878    
## ns(Time, 7) 57.306  6  1.584e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.09$full_model, terms = c("Time")))

em <- emmeans(model.num.09$full_model, ~ Time, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15      -0.0626 0.0853 Inf  -0.734  1.0000
##  0_25       0.0544 0.0732 Inf   0.743  1.0000
##  0_35       0.2319 0.0720 Inf   3.222  0.0102
##  0_45       0.3628 0.0679 Inf   5.345  <.0001
##  0_55       0.3901 0.0658 Inf   5.927  <.0001
##  0_65       0.2943 0.0748 Inf   3.935  0.0007
##  0_75       0.1261 0.0895 Inf   1.409  1.0000
##  0_85       0.0901 0.0950 Inf   0.949  1.0000
## 
## P value adjustment: bonferroni method for 8 tests

The same results as before

FOREARM

# model.num.10 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject),
#                       family = Gamma("identity"), 
#                       data = filter(datas.raw.num, Body == "FOREARM"),
#                       method = "LRT")

Anova(model.num.10$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##               Chisq Df Pr(>Chisq)   
## (Intercept)  0.4698  1   0.493084   
## ns(Time, 7) 19.7598  6   0.003055 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.10$full_model, terms = c("Time")))

em <- emmeans(model.num.10$full_model, ~ Time, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15        0.136 0.0525 Inf   2.597  0.0752
##  0_25        0.182 0.0485 Inf   3.758  0.0014
##  0_35        0.167 0.0509 Inf   3.276  0.0084
##  0_45        0.112 0.0535 Inf   2.094  0.2902
##  0_55        0.143 0.0485 Inf   2.945  0.0259
##  0_65        0.168 0.0507 Inf   3.305  0.0076
##  0_75        0.120 0.0586 Inf   2.053  0.3209
##  0_85        0.194 0.0555 Inf   3.503  0.0037
## 
## P value adjustment: bonferroni method for 8 tests

The same results as before
Exept :
0_55 become significant

ARM

# model.num.11 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject), 
#                       family = Gamma("identity"),
#                       data = filter(datas.raw.num, Body == "ARM"),
#                       method = "LRT")

Anova(model.num.11$full_model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mesure.TP
##               Chisq Df Pr(>Chisq)  
## (Intercept)  0.0054  1    0.94140  
## ns(Time, 7) 16.7440  6    0.01027 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.11$full_model, terms = c("Time")))

em <- emmeans(model.num.11$full_model, ~ Time, 
                            at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")
##  contrast estimate     SE  df z.ratio p.value
##  0_15      -0.0458 0.0732 Inf  -0.626  1.0000
##  0_25       0.0530 0.0618 Inf   0.857  1.0000
##  0_35       0.1219 0.0546 Inf   2.232  0.2049
##  0_45       0.1582 0.0595 Inf   2.661  0.0623
##  0_55       0.1817 0.0571 Inf   3.185  0.0116
##  0_65       0.1383 0.0545 Inf   2.539  0.0889
##  0_75       0.0814 0.0606 Inf   1.345  1.0000
##  0_85       0.0717 0.0681 Inf   1.053  1.0000
## 
## P value adjustment: bonferroni method for 8 tests
# 
# save(model.num.08, model.num.09, model.num.10, model.num.11,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Four_separate.RData")

The same results as before

Power

Reviewers ask effect size or power estimation of our statistics.
This is a tricky question.
We usetwo indirect ways to answer:

Rsquare

Rsquare give an idea of fit quality of our models. So we compare the Rsquare.

library(MuMIn)

body.name <- c("ARM", "FOREARM", "LIP", "CHEEK")
glmer.R2 <- c(r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "ARM")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "FOREARM")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "LIP")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "CHEEK")))[2])

aov.R2<-c(summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "ARM")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "FOREARM")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "LIP")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "CHEEK")))$r.squared)

data.frame(body = body.name,
           glmer = glmer.R2,
           aov = aov.R2) %>%
  mutate(gain = (glmer-aov)/aov)
##      body      glmer         aov     gain
## 1     ARM 0.04769246 0.010683196 3.464250
## 2 FOREARM 0.04216606 0.006509172 5.477945
## 3     LIP 0.05642969 0.020136078 1.802417
## 4   CHEEK 0.05296593 0.016857169 2.142041
data.frame(body = c("LIP & CHEEK", "ARM & FOREARM"),
           glmer = c(r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),data = filter(datas.raw.num, Localization == "Loc1")))[2],
                     r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),data = filter(datas.raw.num, Localization == "Loc2")))[2]),
           aov = c(summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2 ,data = filter(datas.raw.num, Localization == "Loc1")))$r.squared,
                   summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2 ,data = filter(datas.raw.num, Localization == "Loc2"))))$r.squared) %>%
  mutate(gain = (glmer-aov)/aov)
##            body     glmer         aov     gain
## 1   LIP & CHEEK 0.6152456 0.008596999 70.56516
## 2 ARM & FOREARM 0.6094596 0.008596999 69.89213
data.frame(body = c("LIP & CHEEK & ARM & FOREARM"),
           glmer = r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body) + (1|Localization), data = datas.raw.num))[2],
           aov = summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2, data = datas.raw.num))$r.squared) %>%
  mutate(gain = (glmer-aov)/aov)
##                          body     glmer        aov     gain
## 1 LIP & CHEEK & ARM & FOREARM 0.7248271 0.01420603 50.02248

The more we split datas, the less models fit with datas.

Bayes framwork

Another way to check the models quality is to use the same models with bayesian estimation of fit. These model give us some metrics to compare the models quality.

1 model, 4 body parts

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_4.RData")

# model.stan.00 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body) + (1|Localization),                                   family = Gamma("identity"), data = datas.raw.num,
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# 
# save(model.stan.00,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_4.RData")


plot(ggpredict(model.stan.00, terms = c("Time", "Body"))) +
  ylim(c(0.6, 1.25))

emm <- emmeans(model.stan.00, ~ Time | Body, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 80)))


BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.00)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

2 models, 2 body parts each

LIP & CHEEK

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_LIP_CHEEK.RData")

# model.stan.01 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body),                                                       family = Gamma("identity"),
#                             data = filter(datas.raw.num, Localization == "Loc1"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# 
# save(model.stan.01,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_LIP_CHEEK.RData")

plot(ggpredict(model.stan.01, terms = c("Time", "Body")))

emm <- emmeans(model.stan.01, ~ Time | Body, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.01)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

FOREARM & ARM

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_FOREARM_ARM.RData")

# model.stan.02 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body),                                                       family = Gamma("identity"),
#                             data = filter(datas.raw.num, Localization == "Loc2"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# save(model.stan.02,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_FOREARM_ARM.RData")
 
plot(ggpredict(model.stan.02, terms = c("Time", "Body")))

emm <- emmeans(model.stan.02, ~ Time | Body, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.02)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

4 models

LIP

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_LIP.RData")

# model.stan.03 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
#                             data = filter(datas.raw.num, Body == "LIP"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# 
# 
# 
# save(model.stan.03,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_LIP.RData")

plot(ggpredict(model.stan.03, terms = c("Time")))

emm <- emmeans(model.stan.03, ~ Time, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.03)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

CHEEK

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_CHEEK.RData")

# model.stan.04 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
#                             data = filter(datas.raw.num, Body == "CHEEK"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# 
# save(model.stan.04,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_CHEEK.RData")

plot(ggpredict(model.stan.04, terms = c("Time")))

emm <- emmeans(model.stan.04, ~ Time, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.04)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

FOREARM

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_FOREARM.RData")

# model.stan.05 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
#                             data = filter(datas.raw.num, Body == "FOREARM"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# save(model.stan.05,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_FOREARM.RData")

plot(ggpredict(model.stan.05, terms = c("Time")))

emm <- emmeans(model.stan.05, ~ Time, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.05)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

ARM

load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_ARM.RData")

# model.stan.06 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
#                             data = filter(datas.raw.num, Body == "ARM"),
#                             diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# 
# save(model.stan.06,
#      file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_ARM.RData")

plot(ggpredict(model.stan.06, terms = c("Time")))

emm <- emmeans(model.stan.06, ~ Time, 
              at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))

BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.06)

as.data.frame(p_direction(contrast(emm, des_cont))) %>%
  mutate(BF = exp(BF.emm$log_BF)) %>%
  mutate(BayesSign = case_when(BF > 100 ~ "***",
                               BF > 10 & BF < 100 ~ "**",
                               BF > 3 & BF < 10 ~ "*",
                               BF > 1 & BF < 3 ~ ".",
                               BF < 1 ~ "")) -> Tests.Pairs.Results

datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )
rm(Tests.Pairs.Results, emm, BF.emm)

Version

Here, all packages en R version used on this analysis.

print(sessionInfo())
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] MuMIn_1.43.17                bayestestR_0.12.1           
##  [3] bayesplot_1.9.0              brms_2.17.0                 
##  [5] rstanarm_2.21.1              Rcpp_1.0.6                  
##  [7] sjmisc_2.8.9                 sjPlot_2.8.10               
##  [9] glmmTMB_1.1.3                DHARMa_0.4.4                
## [11] ggeffects_1.1.1              ggpubr_0.4.0                
## [13] emmeans_1.7.3                lmerTest_3.1-3              
## [15] multcomp_1.4-16              TH.data_1.0-10              
## [17] MASS_7.3-53                  survival_3.2-7              
## [19] mvtnorm_1.1-1                phia_0.2-1                  
## [21] car_3.0-10                   carData_3.0-4               
## [23] LMERConvenienceFunctions_3.0 RVAideMemoire_0.9-81-2      
## [25] afex_1.1-1                   lme4_1.1-26                 
## [27] Matrix_1.2-18                forcats_0.5.1               
## [29] stringr_1.4.0                dplyr_1.0.8                 
## [31] purrr_0.3.4                  readr_1.4.0                 
## [33] tidyr_1.2.0                  tibble_3.1.6                
## [35] tidyverse_1.3.1              ggplot2_3.3.6               
## [37] openxlsx_4.2.3               readxl_1.3.1                
## [39] readtext_0.80                DT_0.17                     
## [41] epuRate_0.1                  rmarkdown_2.11              
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.1           tidyselect_1.1.2     htmlwidgets_1.5.3   
##   [4] grid_4.0.3           munsell_0.5.0        codetools_0.2-16    
##   [7] effectsize_0.4.4     statmod_1.4.35       miniUI_0.1.1.1      
##  [10] withr_2.4.1          Brobdingnag_1.2-6    colorspace_2.0-0    
##  [13] highr_0.8            knitr_1.39           rstudioapi_0.13     
##  [16] stats4_4.0.3         ggsignif_0.6.1       labeling_0.4.2      
##  [19] LCFdata_2.0          rstan_2.21.1         farver_2.1.0        
##  [22] datawizard_0.4.0     bridgesampling_1.1-2 coda_0.19-4         
##  [25] vctrs_0.4.0          generics_0.1.0       xfun_0.30           
##  [28] R6_2.5.0             markdown_1.1         rmdformats_1.0.2    
##  [31] fields_11.6          gamm4_0.2-6          projpred_2.0.2      
##  [34] logspline_2.1.16     assertthat_0.2.1     promises_1.2.0.1    
##  [37] scales_1.1.1         gtable_0.3.0         processx_3.5.0      
##  [40] spam_2.6-0           sandwich_3.0-0       rlang_1.0.2         
##  [43] Rmisc_1.5            rstatix_0.7.0        TMB_1.7.19          
##  [46] checkmate_2.0.0      broom_0.8.0          inline_0.3.19       
##  [49] yaml_2.2.1           reshape2_1.4.4       abind_1.4-5         
##  [52] modelr_0.1.8         gtsummary_1.6.0      threejs_0.3.3       
##  [55] crosstalk_1.1.1      backports_1.2.1      httpuv_1.5.5        
##  [58] rsconnect_0.8.24     tensorA_0.36.2       tools_4.0.3         
##  [61] bookdown_0.22        ellipsis_0.3.2       RColorBrewer_1.1-2  
##  [64] posterior_1.1.0      jquerylib_0.1.3      ggridges_0.5.3      
##  [67] plyr_1.8.6           base64enc_0.1-3      ps_1.6.0            
##  [70] prettyunits_1.1.1    zoo_1.8-9            haven_2.4.3         
##  [73] fs_1.5.2             magrittr_2.0.3       data.table_1.14.0   
##  [76] colourpicker_1.1.1   reprex_2.0.1         matrixStats_0.58.0  
##  [79] hms_1.0.0            shinyjs_2.0.0        mime_0.10           
##  [82] evaluate_0.15        xtable_1.8-4         shinystan_2.5.0     
##  [85] rio_0.5.26           sjstats_0.18.1       gridExtra_2.3       
##  [88] rstantools_2.2.0     compiler_4.0.3       maps_3.3.0          
##  [91] gt_0.5.0             V8_3.4.2             crayon_1.4.1        
##  [94] minqa_1.2.4          StanHeaders_2.21.0-7 htmltools_0.5.2     
##  [97] mgcv_1.8-33          later_1.1.0.1        RcppParallel_5.1.4  
## [100] lubridate_1.7.10     DBI_1.1.1            sjlabelled_1.1.7    
## [103] dbplyr_2.1.1         broom.helpers_1.7.0  boot_1.3-25         
## [106] cli_3.2.0            parallel_4.0.3       insight_0.17.0      
## [109] dotCall64_1.0-1      igraph_1.2.6         pkgconfig_2.0.3     
## [112] numDeriv_2016.8-1.1  foreign_0.8-80       xml2_1.3.2          
## [115] dygraphs_1.1.1.6     bslib_0.2.4          estimability_1.3    
## [118] rvest_1.0.0          distributional_0.2.2 callr_3.7.0         
## [121] digest_0.6.27        parameters_0.14.0    cellranger_1.1.0    
## [124] curl_4.3             commonmark_1.7       shiny_1.6.0         
## [127] gtools_3.8.2         nloptr_1.2.2.2       lifecycle_1.0.1     
## [130] nlme_3.1-149         jsonlite_1.7.2       fansi_0.4.2         
## [133] pillar_1.7.0         lattice_0.20-41      loo_2.4.1           
## [136] fastmap_1.1.0        httr_1.4.2           pkgbuild_1.2.0      
## [139] glue_1.6.2           xts_0.12.1           zip_2.1.1           
## [142] shinythemes_1.2.0    stringi_1.5.3        sass_0.4.1          
## [145] performance_0.9.0